knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
fig.height = 4, fig.width = 7)
library(terra)
library(sf)
library(oharac)
library(data.table)
library(tidyverse)
library(here)
library(cowplot)
source(here('common_fxns.R'))Sensitivity analysis to determine how drastically scores change with different values of functional traits.
In script 4a, we calculated sensitivity around functional diversity metrics and per-stressor impacts. In 4b, we aggregate results across all stressors for each randomized iteration, to determine the distribution of mean values for cumulative impacts. Here, we map them out and make cool figures.
For each randomized trait, plot the variation in functional diversity metrics for the sampled points, with a nearest-neighbor interpolation to help spot regional variation. We will examine several stats: absolute and relative difference (diff / mean), 95% confidence interval, coefficient of variance (std dev divided by mean)
ocean_r <- rast(here('_spatial/ocean_area_mol.tif'))
cell_xy <- ocean_r %>%
setValues(1:ncell(.)) %>%
setNames('cell_id') %>%
as.data.frame(xy = TRUE)fd_fs <- list.files(here_anx('sens_analysis/funct_div_summary'),
pattern = '.csv', full.names = TRUE)
fd_ts <- basename(fd_fs) %>% str_remove_all('.+_shuffle_|.csv')
fd_df <- lapply(fd_fs, fread) %>%
setNames(fd_ts) %>%
rbindlist(idcol = 't') %>%
mutate(n_fe_diff = n_fe_mean - n_fe,
fv_diff = fv_mean - fv,
n_fe_ci = (n_fe_95hi - n_fe_95lo),
fv_ci = (fv_95hi - fv_95lo),
n_fe_cv = n_fe_sdev / n_fe_mean,
fv_cv = fv_sdev / fv_mean) %>%
rename(n_fe_value = n_fe, fv_value = fv) %>%
pivot_longer(cols = -c(cell_id, t, nspp),
names_to = 'var', values_to = 'val') %>%
mutate(param = str_extract(var, 'n_fe|fv'),
stat = str_extract(var, 'mean|sdev|ci|cv|diff|95..|value')) %>%
select(-var) %>%
left_join(cell_xy, by = 'cell_id')How does functional vulnerability vary as traits are randomized?
land_sf <- read_sf(here('_spatial/ne_10m_land/ne_10m_land_no_casp.shp')) %>%
st_transform(st_crs(ocean_r))
globe_bbox <- rbind(c(-180, -90), c(-180, 90),
c(180, 90), c(180, -90), c(-180, -90))
globe_border <- st_polygon(list(globe_bbox)) %>%
st_sfc(crs = 4326) %>%
st_sf(data.frame(rgn = 'globe', geom = .)) %>%
smoothr::densify(max_distance = 0.5) %>%
st_transform(crs = crs(ocean_r))map_param <- function(df, s, p, r = ocean_r, sig = FALSE, thresh = NULL) {
val_mtx <- df %>%
filter(stat == s & param == p) %>%
select(x, y, z = val) %>%
as.matrix()
val_rast <- interpNear(r, val_mtx, radius = 500000) %>%
mask(r)
if(sig) {
sig_mtx <- df %>%
filter(param == p & stat %in% c('95hi', '95lo', 'value')) %>%
pivot_wider(names_from = stat, values_from = val) %>%
mutate(sig = ifelse(between(value, `95lo`, `95hi`), 0, 1)) %>%
select(x, y, z = sig) %>%
as.matrix()
sig_rast <- interpNear(r, sig_mtx, radius = 50000) %>%
mask(r)
val_rast[sig_rast == 0] <- NA
}
if(!is.null(thresh)) val_rast[abs(val_rast) < thresh] <- NA
return(val_rast)
}
finish_map <- function(p, div = TRUE) {
### div is whether the palette should be a diverging palette (TRUE) or
### sequential (FALSE)
p <- p +
### continents:
geom_sf(data = land_sf,
fill = 'grey96', color = 'grey40',
size = .10) +
geom_sf(data = globe_border,
fill = NA, color = 'grey70',
size = .1) +
scale_x_continuous(expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
facet_wrap(~ trait) +
theme_void()
if(div) {
### diverging palette
p <- p +
scale_fill_gradient2(low = '#d01c8b', mid = '#ffffdf', high = '#4dac26')
} else {
### sequential palette
p <- p + scale_fill_viridis_c()
}
return(p)
}fv_fig_f <- '_figs/fv_diff_by_trait_map.png'
if(!file.exists(fv_fig_f)) {
fv_map_list <- vector('list', length = length(fd_ts)) %>%
setNames(fd_ts)
for(trait in fd_ts) {
# trait <- fd_ts[1]
fv_map_list[[trait]] <- map_param(df = fd_df %>% filter(t == trait),
p = 'fv', s = 'diff')
}
fv_map_stack <- rast(fv_map_list)
mean_sd_vals <- as.data.frame(fv_map_stack) %>%
summarize(across(everything(), .fns = list(mean = mean, sd = sd), .names = '{.col}_{.fn}'))
# len_mean len_sd mob_mean mob_sd trp_mean trp_sd wcol_mean wcol_sd
# 1 0.03633012 0.04194202 0.03743112 0.04292291 0.05350717 0.04770103 0.02206068 0.02724206
fv_map_df <- as.data.frame(fv_map_stack, xy = TRUE) %>%
pivot_longer(-c(x, y), names_to = 'trait', values_to = 'delta_fv') %>%
mutate(trait = case_when(trait == 'len' ~ 'body length',
trait == 'mob' ~ 'adult mobility',
trait == 'trp' ~ 'trophic level',
trait == 'wcol' ~ 'water column position'))
p <- ggplot() +
geom_raster(data = fv_map_df, aes(x, y, fill = delta_fv)) +
labs(fill = 'ΔFV')
p <- finish_map(p)
ggsave(fv_fig_f, width = 6.5, height = 3, dpi = 300)
}
knitr::include_graphics(fv_fig_f)fv2_fig_f <- '_figs/fv_cv_by_trait_map.png'
if(!file.exists(fv2_fig_f)) {
fv_map_list <- vector('list', length = length(fd_ts)) %>%
setNames(fd_ts)
for(trait in fd_ts) {
fv_map_list[[trait]] <- map_param(df = fd_df %>% filter(t == trait),
p = 'fv', s = 'cv')
}
fv_map_stack <- rast(fv_map_list)
fv_map_df <- as.data.frame(fv_map_stack, xy = TRUE) %>%
pivot_longer(-c(x, y), names_to = 'trait', values_to = 'cv_fv') %>%
mutate(trait = case_when(trait == 'len' ~ 'body length',
trait == 'mob' ~ 'adult mobility',
trait == 'trp' ~ 'trophic level',
trait == 'wcol' ~ 'water column position'))
p <- ggplot() +
geom_raster(data = fv_map_df, aes(x, y, fill = cv_fv)) +
labs(fill = 'CV(FV)')
p <- finish_map(p, div = FALSE)
ggsave(fv2_fig_f, width = 6.5, height = 3, dpi = 300)
}
knitr::include_graphics(fv2_fig_f)How do functional diversity metrics (FV and number of FEs) vary as traits are randomized, relative to the number of species in a given pixel?
plot_f <- '_figs/fv_stat_by_trait_scatter.png'
if(!file.exists(plot_f)) {
fv_plot_df <- fd_df %>%
filter(stat %in% c('mean', 'cv', 'diff')) %>%
filter(param == 'fv') %>%
mutate(metric = paste(stat, param)) %>%
mutate(metric = str_replace(metric, 'fv', '(FV)'),
metric = str_replace(metric, 'mean ', 'mean'),
metric = str_replace(metric, 'cv ', 'CV'),
metric = str_replace(metric, 'diff ', 'Δ')) %>%
mutate(t = case_when(t == 'len' ~ 'body length',
t == 'mob' ~ 'adult mobility',
t == 'trp' ~ 'trophic level',
t == 'wcol' ~ 'water column position')) %>%
sample_frac(.1)
p <- ggplot(fv_plot_df, aes(x = nspp, y = val)) +
geom_point(alpha = .2) +
theme_ohara() +
facet_grid(metric ~ t, scales = 'free_y') +
labs(x = 'Species richness') +
theme(axis.title.y = element_blank())
ggsave(plot_f, width = 6, height = 3, dpi = 300)
}
knitr::include_graphics(plot_f)plot_f <- '_figs/n_fe_stat_by_trait_scatter.png'
if(!file.exists(plot_f)) {
fe_plot_df <- fd_df %>%
filter(stat %in% c('mean', 'cv', 'diff')) %>%
filter(param == 'n_fe') %>%
mutate(metric = paste(stat, param)) %>%
mutate(metric = str_replace(metric, 'n_fe', '(# FE)'),
metric = str_replace(metric, 'mean ', 'mean'),
metric = str_replace(metric, 'cv ', 'CV'),
metric = str_replace(metric, 'diff ', 'Δ')) %>%
mutate(t = case_when(t == 'len' ~ 'body length',
t == 'mob' ~ 'adult mobility',
t == 'trp' ~ 'trophic level',
t == 'wcol' ~ 'water column position')) %>%
sample_frac(.1)
p <- ggplot(fe_plot_df, aes(x = nspp, y = val)) +
geom_point(alpha = .2) +
theme_ohara() +
facet_grid(metric ~ t, scales = 'free_y') +
labs(x = 'Species richness') +
theme(axis.title.y = element_blank())
ggsave(plot_f, width = 6, height = 3, dpi = 300)
}
knitr::include_graphics(plot_f)How does cumulative impact vary as traits are randomized?
chi_fs <- list.files(here_anx('sens_analysis/impact_summary'),
pattern = 'chi_.+.csv', full.names = TRUE)
chi_ts <- basename(chi_fs) %>% str_remove_all('.+_shuffle_|.csv')
chi_true_df <- rast(here('_output/cumulative_impact_maps/chi_funct_entity.tif')) %>%
as.data.frame(xy = TRUE) %>%
oharac::dt_join(cell_xy, by = c('x', 'y'), type = 'left') %>%
filter(cell_id %in% fd_df$cell_id)
nspp_df <- fd_df %>%
select(cell_id, nspp) %>%
distinct()
chi_df <- lapply(chi_fs, fread) %>%
setNames(chi_ts) %>%
rbindlist(idcol = 't') %>%
left_join(chi_true_df, by = 'cell_id') %>%
left_join(nspp_df, by = 'cell_id') %>%
mutate(chi_diff = chi_mean - chi_fe,
chi_ci = (chi_95hi - chi_95lo),
chi_cv = chi_sdev / chi_mean) %>%
select(cell_id, t, nspp,
chi_mean, chi_diff,
chi_ci, chi_cv) %>%
pivot_longer(cols = -c(cell_id, t, nspp),
names_to = 'var', values_to = 'val') %>%
mutate(param = 'chi',
stat = str_extract(var, 'mean|sdev|ci|cv|diff')) %>%
select(-var) %>%
left_join(cell_xy, by = 'cell_id')chi_fig_f <- '_figs/chi_diff_by_trait_map.png'
if(!file.exists(chi_fig_f)) {
chi_map_list <- vector('list', length = length(chi_ts)) %>%
setNames(chi_ts)
for(trait in fd_ts) {
chi_map_list[[trait]] <- map_param(df = chi_df %>% filter(t == trait),
p = 'chi', s = 'diff')
}
chi_map_stack <- rast(chi_map_list)
mean_sd_vals <- as.data.frame(chi_map_stack) %>%
summarize(across(everything(), .fns = list(mean = mean, sd = sd), .names = '{.col}_{.fn}'))
# len_mean len_sd mob_mean mob_sd trp_mean trp_sd wcol_mean wcol_sd
# 1 -0.00801914 0.01990471 0.0008610911 0.01656777 0.003189031 0.01596556 -0.006297416 0.02280766
chi_map_df <- as.data.frame(chi_map_stack, xy = TRUE) %>%
pivot_longer(-c(x, y), names_to = 'trait', values_to = 'delta_chi') %>%
mutate(trait = case_when(trait == 'len' ~ 'body length',
trait == 'mob' ~ 'adult mobility',
trait == 'trp' ~ 'trophic level',
trait == 'wcol' ~ 'water column position'))
p <- ggplot() +
geom_raster(data = chi_map_df, aes(x, y, fill = delta_chi)) +
labs(fill = 'ΔCHI')
p <- finish_map(p)
ggsave(chi_fig_f, width = 6.5, height = 3, dpi = 300)
}
knitr::include_graphics(chi_fig_f)How do CHI stats vary as traits are randomized, relative to the number of species in a given pixel?
plot_f <- '_figs/chi_stat_by_trait_scatter.png'
if(!file.exists(plot_f)) {
chi_plot_df <- chi_df %>%
filter(stat %in% c('mean', 'cv', 'diff')) %>%
mutate(metric = paste(stat, param)) %>%
mutate(metric = str_replace(metric, 'chi', '(CHI)'),
metric = str_replace(metric, 'mean ', 'mean'),
metric = str_replace(metric, 'cv ', 'CV'),
metric = str_replace(metric, 'diff ', 'Δ')) %>%
mutate(t = case_when(t == 'len' ~ 'body length',
t == 'mob' ~ 'adult mobility',
t == 'trp' ~ 'trophic level',
t == 'wcol' ~ 'water column position')) %>%
sample_frac(.1)
p <- ggplot(chi_plot_df, aes(x = nspp, y = val)) +
geom_point(alpha = .2) +
theme_ohara() +
facet_grid(metric ~ t, scales = 'free_y') +
labs(x = 'Species richness') +
theme(axis.title.y = element_blank())
ggsave(plot_f, width = 6, height = 3, dpi = 300)
}
knitr::include_graphics(plot_f)Combine the FV variation and CHI variation into a large panel.
p <- ggdraw() +
draw_image(fv_fig_f, x = 0.0, y = 0.51, width = 1.0, height = 0.49) +
draw_image(chi_fig_f, x = 0.0, y = 0.0, width = 1.0, height = 0.49) +
draw_label('A.', x = .07, y = .995, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('B.', x = 0.5, y = .995, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('C.', x = .07, y = .748, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('D.', x = 0.5, y = .748, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('E.', x = .07, y = .485, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('F.', x = 0.5, y = .485, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('G.', x = .07, y = .235, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold') +
draw_label('H.', x = 0.5, y = .235, hjust = 0, vjust = 1,
size = 7, color = 'black', fontfamily = 'Helvetica', fontface = 'bold')
si_fig_file <- here('5_ms_figs/figS6_sens_analysis.png')
ggsave(si_fig_file, height = 11.5, width = 12, units = 'cm', dpi = 300)
knitr::include_graphics(si_fig_file)